home *** CD-ROM | disk | FTP | other *** search
/ Sprite 1984 - 1993 / Sprite 1984 - 1993.iso / src / lib / m / support.c < prev    next >
C/C++ Source or Header  |  1990-02-16  |  17KB  |  507 lines
  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms are permitted
  6.  * provided that this notice is preserved and that due credit is given
  7.  * to the University of California at Berkeley. The name of the University
  8.  * may not be used to endorse or promote products derived from this
  9.  * software without specific prior written permission. This software
  10.  * is provided ``as is'' without express or implied warranty.
  11.  *
  12.  * All recipients should regard themselves as participants in an ongoing
  13.  * research project and hence should feel obligated to report their
  14.  * experiences (good or bad) with these elementary function codes, using
  15.  * the sendbug(8) program, to the authors.
  16.  */
  17.  
  18. #ifndef lint
  19. static char sccsid[] = "@(#)support.c    5.2 (Berkeley) 4/29/88";
  20. #endif /* not lint */
  21.  
  22. /* 
  23.  * Some IEEE standard 754 recommended functions and remainder and sqrt for 
  24.  * supporting the C elementary functions.
  25.  ******************************************************************************
  26.  * WARNING:
  27.  *      These codes are developed (in double) to support the C elementary
  28.  * functions temporarily. They are not universal, and some of them are very
  29.  * slow (in particular, drem and sqrt is extremely inefficient). Each 
  30.  * computer system should have its implementation of these functions using 
  31.  * its own assembler.
  32.  ******************************************************************************
  33.  *
  34.  * IEEE 754 required operations:
  35.  *     drem(x,p) 
  36.  *              returns  x REM y  =  x - [x/y]*y , where [x/y] is the integer
  37.  *              nearest x/y; in half way case, choose the even one.
  38.  *     sqrt(x) 
  39.  *              returns the square root of x correctly rounded according to 
  40.  *        the rounding mod.
  41.  *
  42.  * IEEE 754 recommended functions:
  43.  * (a) copysign(x,y) 
  44.  *              returns x with the sign of y. 
  45.  * (b) scalb(x,N) 
  46.  *              returns  x * (2**N), for integer values N.
  47.  * (c) logb(x) 
  48.  *              returns the unbiased exponent of x, a signed integer in 
  49.  *              double precision, except that logb(0) is -INF, logb(INF) 
  50.  *              is +INF, and logb(NAN) is that NAN.
  51.  * (d) finite(x) 
  52.  *              returns the value TRUE if -INF < x < +INF and returns 
  53.  *              FALSE otherwise.
  54.  *
  55.  *
  56.  * CODED IN C BY K.C. NG, 11/25/84;
  57.  * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
  58.  */
  59.  
  60.  
  61. #if defined(vax)||defined(tahoe)      /* VAX D format */
  62. #include <errno.h>
  63.     static unsigned short msign=0x7fff , mexp =0x7f80 ;
  64.     static short  prep1=57, gap=7, bias=129           ;   
  65.     static double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
  66. #else    /* defined(vax)||defined(tahoe) */
  67.     static unsigned short msign=0x7fff, mexp =0x7ff0  ;
  68.     static short prep1=54, gap=4, bias=1023           ;
  69.     static double novf=1.7E308, nunf=3.0E-308,zero=0.0;
  70. #endif    /* defined(vax)||defined(tahoe) */
  71.  
  72. double scalb(x,N)
  73. double x; int N;
  74. {
  75.         int k;
  76.         double scalb();
  77.  
  78. #if defined(national) || defined (ds3100)
  79.         unsigned short *px=(unsigned short *) &x + 3;
  80. #else    /* national  || ds3100 */
  81.         unsigned short *px=(unsigned short *) &x;
  82. #endif    /* national */
  83.  
  84.         if( x == zero )  return(x); 
  85.  
  86. #if defined(vax)||defined(tahoe)
  87.         if( (k= *px & mexp ) != ~msign ) {
  88.             if (N < -260)
  89.         return(nunf*nunf);
  90.         else if (N > 260) {
  91.         extern double infnan(),copysign();
  92.         return(copysign(infnan(ERANGE),x));
  93.         }
  94. #else    /* defined(vax)||defined(tahoe) */
  95.         if( (k= *px & mexp ) != mexp ) {
  96.             if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf);
  97.             if( k == 0 ) {
  98.                  x *= scalb(1.0,(int)prep1);  N -= prep1; return(scalb(x,N));}
  99. #endif    /* defined(vax)||defined(tahoe) */
  100.  
  101.             if((k = (k>>gap)+ N) > 0 )
  102.                 if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
  103.                 else x=novf+novf;               /* overflow */
  104.             else
  105.                 if( k > -prep1 ) 
  106.                                         /* gradual underflow */
  107.                     {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);}
  108.                 else
  109.                 return(nunf*nunf);
  110.             }
  111.         return(x);
  112. }
  113.  
  114.  
  115. double copysign(x,y)
  116. double x,y;
  117. {
  118. #ifdef national
  119.         unsigned short  *px=(unsigned short *) &x+3,
  120.                         *py=(unsigned short *) &y+3;
  121. #else    /* national */
  122.         unsigned short  *px=(unsigned short *) &x,
  123.                         *py=(unsigned short *) &y;
  124. #endif    /* national */
  125.  
  126. #if defined(vax)||defined(tahoe)
  127.         if ( (*px & mexp) == 0 ) return(x);
  128. #endif    /* defined(vax)||defined(tahoe) */
  129.  
  130.         *px = ( *px & msign ) | ( *py & ~msign );
  131.         return(x);
  132. }
  133.  
  134. double logb(x)
  135. double x; 
  136. {
  137.  
  138. #ifdef national
  139.         short *px=(short *) &x+3, k;
  140. #else    /* national */
  141.         short *px=(short *) &x, k;
  142. #endif    /* national */
  143.  
  144. #if defined(vax)||defined(tahoe)
  145.         return (int)(((*px&mexp)>>gap)-bias);
  146. #else    /* defined(vax)||defined(tahoe) */
  147.         if( (k= *px & mexp ) != mexp )
  148.             if ( k != 0 )
  149.                 return ( (k>>gap) - bias );
  150.             else if( x != zero)
  151.                 return ( -1022.0 );
  152.             else        
  153.                 return(-(1.0/zero));    
  154.         else if(x != x)
  155.             return(x);
  156.         else
  157.             {*px &= msign; return(x);}
  158. #endif    /* defined(vax)||defined(tahoe) */
  159. }
  160.  
  161. finite(x)
  162. double x;    
  163. {
  164. #if defined(vax)||defined(tahoe)
  165.         return(1);
  166. #else    /* defined(vax)||defined(tahoe) */
  167. #ifdef national
  168.         return( (*((short *) &x+3 ) & mexp ) != mexp );
  169. #else    /* national */
  170.         return( (*((short *) &x ) & mexp ) != mexp );
  171. #endif    /* national */
  172. #endif    /* defined(vax)||defined(tahoe) */
  173. }
  174.  
  175. double drem(x,p)
  176. double x,p;
  177. {
  178.         short sign;
  179.         double hp,dp,tmp,drem(),scalb();
  180.         unsigned short  k; 
  181. #ifdef national
  182.         unsigned short
  183.               *px=(unsigned short *) &x  +3, 
  184.               *pp=(unsigned short *) &p  +3,
  185.               *pd=(unsigned short *) &dp +3,
  186.               *pt=(unsigned short *) &tmp+3;
  187. #else    /* national */
  188.         unsigned short
  189.               *px=(unsigned short *) &x  , 
  190.               *pp=(unsigned short *) &p  ,
  191.               *pd=(unsigned short *) &dp ,
  192.               *pt=(unsigned short *) &tmp;
  193. #endif    /* national */
  194.  
  195.         *pp &= msign ;
  196.  
  197. #if defined(vax)||defined(tahoe)
  198.         if( ( *px & mexp ) == ~msign )    /* is x a reserved operand? */
  199. #else    /* defined(vax)||defined(tahoe) */
  200.         if( ( *px & mexp ) == mexp )
  201. #endif    /* defined(vax)||defined(tahoe) */
  202.         return  (x-p)-(x-p);    /* create nan if x is inf */
  203.     if (p == zero) {
  204. #if defined(vax)||defined(tahoe)
  205.         extern double infnan();
  206.         return(infnan(EDOM));
  207. #else    /* defined(vax)||defined(tahoe) */
  208.         return zero/zero;
  209. #endif    /* defined(vax)||defined(tahoe) */
  210.     }
  211.  
  212. #if defined(vax)||defined(tahoe)
  213.         if( ( *pp & mexp ) == ~msign )    /* is p a reserved operand? */
  214. #else    /* defined(vax)||defined(tahoe) */
  215.         if( ( *pp & mexp ) == mexp )
  216. #endif    /* defined(vax)||defined(tahoe) */
  217.         { if (p != p) return p; else return x;}
  218.  
  219.         else  if ( ((*pp & mexp)>>gap) <= 1 ) 
  220.                 /* subnormal p, or almost subnormal p */
  221.             { double b; b=scalb(1.0,(int)prep1);
  222.               p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);}
  223.         else  if ( p >= novf/2)
  224.             { p /= 2 ; x /= 2; return(drem(x,p)*2);}
  225.         else 
  226.             {
  227.                 dp=p+p; hp=p/2;
  228.                 sign= *px & ~msign ;
  229.                 *px &= msign       ;
  230.                 while ( x > dp )
  231.                     {
  232.                         k=(*px & mexp) - (*pd & mexp) ;
  233.                         tmp = dp ;
  234.                         *pt += k ;
  235.  
  236. #if defined(vax)||defined(tahoe)
  237.                         if( x < tmp ) *pt -= 128 ;
  238. #else    /* defined(vax)||defined(tahoe) */
  239.                         if( x < tmp ) *pt -= 16 ;
  240. #endif    /* defined(vax)||defined(tahoe) */
  241.  
  242.                         x -= tmp ;
  243.                     }
  244.                 if ( x > hp )
  245.                     { x -= p ;  if ( x >= hp ) x -= p ; }
  246.  
  247. #if defined(vax)||defined(tahoe)
  248.         if (x)
  249. #endif    /* defined(vax)||defined(tahoe) */
  250.             *px ^= sign;
  251.                 return( x);
  252.  
  253.             }
  254. }
  255. double sqrt(x)
  256. double x;
  257. {
  258.         double q,s,b,r;
  259.         double logb(),scalb();
  260.         double t,zero=0.0;
  261.         int m,n,i,finite();
  262. #if defined(vax)||defined(tahoe)
  263.         int k=54;
  264. #else    /* defined(vax)||defined(tahoe) */
  265.         int k=51;
  266. #endif    /* defined(vax)||defined(tahoe) */
  267.  
  268.     /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
  269.         if(x!=x||x==zero) return(x);
  270.  
  271.     /* sqrt(negative) is invalid */
  272.         if(x<zero) {
  273. #if defined(vax)||defined(tahoe)
  274.         extern double infnan();
  275.         return (infnan(EDOM));    /* NaN */
  276. #else    /* defined(vax)||defined(tahoe) */
  277.         return(zero/zero);
  278. #endif    /* defined(vax)||defined(tahoe) */
  279.     }
  280.  
  281.     /* sqrt(INF) is INF */
  282.         if(!finite(x)) return(x);               
  283.  
  284.     /* scale x to [1,4) */
  285.         n=logb(x);
  286.         x=scalb(x,-n);
  287.         if((m=logb(x))!=0) x=scalb(x,-m);       /* subnormal number */
  288.         m += n; 
  289.         n = m/2;
  290.         if((n+n)!=m) {x *= 2; m -=1; n=m/2;}
  291.  
  292.     /* generate sqrt(x) bit by bit (accumulating in q) */
  293.             q=1.0; s=4.0; x -= 1.0; r=1;
  294.             for(i=1;i<=k;i++) {
  295.                 t=s+1; x *= 4; r /= 2;
  296.                 if(t<=x) {
  297.                     s=t+t+2, x -= t; q += r;}
  298.                 else
  299.                     s *= 2;
  300.                 }
  301.             
  302.     /* generate the last bit and determine the final rounding */
  303.             r/=2; x *= 4; 
  304.             if(x==zero) goto end; 100+r; /* trigger inexact flag */
  305.             if(s<x) {
  306.                 q+=r; x -=s; s += 2; s *= 2; x *= 4;
  307.                 t = (x-s)-5; 
  308.                 b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */
  309.                 b=1.0+r/4;   if(b>1.0) t=1;    /* b>1 : Round-to-(+INF) */
  310.                 if(t>=0) q+=r; }          /* else: Round-to-nearest */
  311.             else { 
  312.                 s *= 2; x *= 4; 
  313.                 t = (x-s)-1; 
  314.                 b=1.0+3*r/4; if(b==1.0) goto end;
  315.                 b=1.0+r/4;   if(b>1.0) t=1;
  316.                 if(t>=0) q+=r; }
  317.             
  318. end:        return(scalb(q,n));
  319. }
  320.  
  321. #if 0
  322. /* DREM(X,Y)
  323.  * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
  324.  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
  325.  * INTENDED FOR ASSEMBLY LANGUAGE
  326.  * CODED IN C BY K.C. NG, 3/23/85, 4/8/85.
  327.  *
  328.  * Warning: this code should not get compiled in unless ALL of
  329.  * the following machine-dependent routines are supplied.
  330.  * 
  331.  * Required machine dependent functions (not on a VAX):
  332.  *     swapINX(i): save inexact flag and reset it to "i"
  333.  *     swapENI(e): save inexact enable and reset it to "e"
  334.  */
  335.  
  336. double drem(x,y)    
  337. double x,y;
  338. {
  339.  
  340. #ifdef national        /* order of words in floating point number */
  341.     static n0=3,n1=2,n2=1,n3=0;
  342. #else /* VAX, SUN, ZILOG, TAHOE */
  343.     static n0=0,n1=1,n2=2,n3=3;
  344. #endif
  345.  
  346.         static unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390;
  347.     static double zero=0.0;
  348.     double hy,y1,t,t1;
  349.     short k;
  350.     long n;
  351.     int i,e; 
  352.     unsigned short xexp,yexp, *px  =(unsigned short *) &x  , 
  353.                   nx,nf,      *py  =(unsigned short *) &y  ,
  354.                   sign,      *pt  =(unsigned short *) &t  ,
  355.                         *pt1 =(unsigned short *) &t1 ;
  356.  
  357.     xexp = px[n0] & mexp ;    /* exponent of x */
  358.     yexp = py[n0] & mexp ;    /* exponent of y */
  359.     sign = px[n0] &0x8000;    /* sign of x     */
  360.  
  361. /* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */
  362.     if(x!=x) return(x); if(y!=y) return(y);         /* x or y is NaN */
  363.     if( xexp == mexp )   return(zero/zero);      /* x is INF */
  364.     if(y==zero) return(y/y);
  365.  
  366. /* save the inexact flag and inexact enable in i and e respectively
  367.  * and reset them to zero
  368.  */
  369.     i=swapINX(0);    e=swapENI(0);    
  370.  
  371. /* subnormal number */
  372.     nx=0;
  373.     if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;}
  374.  
  375. /* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */
  376.     if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;}
  377.  
  378.     nf=nx;
  379.     py[n0] &= 0x7fff;    
  380.     px[n0] &= 0x7fff;
  381.  
  382. /* mask off the least significant 27 bits of y */
  383.     t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t;
  384.  
  385. /* LOOP: argument reduction on x whenever x > y */
  386. loop:
  387.     while ( x > y )
  388.     {
  389.         t=y;
  390.         t1=y1;
  391.         xexp=px[n0]&mexp;      /* exponent of x */
  392.         k=xexp-yexp-m25;
  393.         if(k>0)     /* if x/y >= 2**26, scale up y so that x/y < 2**26 */
  394.         {pt[n0]+=k;pt1[n0]+=k;}
  395.         n=x/t; x=(x-n*t1)-n*(t-t1);
  396.     }    
  397.     /* end while (x > y) */
  398.  
  399.     if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;}
  400.  
  401. /* final adjustment */
  402.  
  403.     hy=y/2.0;
  404.     if(x>hy||((x==hy)&&n%2==1)) x-=y; 
  405.     px[n0] ^= sign;
  406.     if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;}
  407.  
  408. /* restore inexact flag and inexact enable */
  409.     swapINX(i); swapENI(e);    
  410.  
  411.     return(x);    
  412. }
  413. #endif
  414.  
  415. #if 0
  416. /* SQRT
  417.  * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT
  418.  * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE
  419.  * CODED IN C BY K.C. NG, 3/22/85.
  420.  *
  421.  * Warning: this code should not get compiled in unless ALL of
  422.  * the following machine-dependent routines are supplied.
  423.  * 
  424.  * Required machine dependent functions:
  425.  *     swapINX(i)  ...return the status of INEXACT flag and reset it to "i"
  426.  *     swapRM(r)   ...return the current Rounding Mode and reset it to "r"
  427.  *     swapENI(e)  ...return the status of inexact enable and reset it to "e"
  428.  *     addc(t)     ...perform t=t+1 regarding t as a 64 bit unsigned integer
  429.  *     subc(t)     ...perform t=t-1 regarding t as a 64 bit unsigned integer
  430.  */
  431.  
  432. static unsigned long table[] = {
  433. 0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740,
  434. 58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478,
  435. 21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, };
  436.  
  437. double newsqrt(x)
  438. double x;
  439. {
  440.         double y,z,t,addc(),subc(),b54=134217728.*134217728.; /* b54=2**54 */
  441.         long mx,scalx,mexp=0x7ff00000;
  442.         int i,j,r,e,swapINX(),swapRM(),swapENI();       
  443.         unsigned long *py=(unsigned long *) &y   ,
  444.                       *pt=(unsigned long *) &t   ,
  445.                       *px=(unsigned long *) &x   ;
  446. #ifdef national         /* ordering of word in a floating point number */
  447.         int n0=1, n1=0; 
  448. #else
  449.         int n0=0, n1=1; 
  450. #endif
  451. /* Rounding Mode:  RN ...round-to-nearest 
  452.  *                 RZ ...round-towards 0
  453.  *                 RP ...round-towards +INF
  454.  *           RM ...round-towards -INF
  455.  */
  456.         int RN=0,RZ=1,RP=2,RM=3;/* machine dependent: work on a Zilog Z8070
  457.                                  * and a National 32081 & 16081
  458.                                  */
  459.  
  460. /* exceptions */
  461.     if(x!=x||x==0.0) return(x);  /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
  462.     if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */
  463.         if((mx=px[n0]&mexp)==mexp) return(x);  /* sqrt(+INF) is +INF */
  464.  
  465. /* save, reset, initialize */
  466.         e=swapENI(0);   /* ...save and reset the inexact enable */
  467.         i=swapINX(0);   /* ...save INEXACT flag */
  468.         r=swapRM(RN);   /* ...save and reset the Rounding Mode to RN */
  469.         scalx=0;
  470.  
  471. /* subnormal number, scale up x to x*2**54 */
  472.         if(mx==0) {x *= b54 ; scalx-=0x01b00000;}
  473.  
  474. /* scale x to avoid intermediate over/underflow:
  475.  * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */
  476.         if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;}
  477.         if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;}
  478.  
  479. /* magic initial approximation to almost 8 sig. bits */
  480.         py[n0]=(px[n0]>>1)+0x1ff80000;
  481.         py[n0]=py[n0]-table[(py[n0]>>15)&31];
  482.  
  483. /* Heron's rule once with correction to improve y to almost 18 sig. bits */
  484.         t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0;
  485.  
  486. /* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */
  487.         t=y*y; z=t;  pt[n0]+=0x00100000; t+=z; z=(x-z)*y; 
  488.         t=z/(t+x) ;  pt[n0]+=0x00100000; y+=t;
  489.  
  490. /* twiddle last bit to force y correctly rounded */ 
  491.         swapRM(RZ);     /* ...set Rounding Mode to round-toward-zero */
  492.         swapINX(0);     /* ...clear INEXACT flag */
  493.         swapENI(e);     /* ...restore inexact enable status */
  494.         t=x/y;          /* ...chopped quotient, possibly inexact */
  495.         j=swapINX(i);   /* ...read and restore inexact flag */
  496.         if(j==0) { if(t==y) goto end; else t=subc(t); }  /* ...t=t-ulp */
  497.         b54+0.1;        /* ..trigger inexact flag, sqrt(x) is inexact */
  498.         if(r==RN) t=addc(t);            /* ...t=t+ulp */
  499.         else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */
  500.         y=y+t;                          /* ...chopped sum */
  501.         py[n0]=py[n0]-0x00100000;       /* ...correctly rounded sqrt(x) */
  502. end:    py[n0]=py[n0]+scalx;            /* ...scale back y */
  503.         swapRM(r);                      /* ...restore Rounding Mode */
  504.         return(y);
  505. }
  506. #endif
  507.